library(fpp3)
── Attaching packages ────────────────────────────────────────────────────── fpp3 0.1 ──
✔ tibble      2.1.1     ✔ tsibble     0.8.5
✔ dplyr       0.8.3     ✔ tsibbledata 0.1.0
✔ tidyr       1.0.0     ✔ feasts      0.1.2
✔ lubridate   1.7.4     ✔ fable       0.1.1
✔ ggplot2     3.1.1     
package ‘fabletools’ was built under R version 3.6.2── Conflicts ───────────────────────────────────────────────────────── fpp3_conflicts ──
✖ lubridate::date()       masks base::date()
✖ dplyr::filter()         masks stats::filter()
✖ tsibble::id()           masks dplyr::id()
✖ tsibble::interval()     masks lubridate::interval()
✖ dplyr::lag()            masks stats::lag()
✖ tsibble::new_interval() masks lubridate::new_interval()

Exponential smoothing

see pdf notes for details on additive and multiplicative models.

ETS(error, trend, season)

Each error, trend and season can be additive/multiplicative/none, or damped.

report(fit)
Series: Pop 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.3266366 

  Initial states:
        l         b
 10.05414 0.2224818

  sigma^2:  0.0041

      AIC      AICc       BIC 
-76.98569 -75.83184 -66.68347 

ETS(A,A,N) : A stands for “additive”, M for “multiplicative”, and N for “none”. Damped is “Ad” or “Md” for additive and multiplicative respectively.

\(\sigma^2\) is the variance of the residuals

Default chooses best (which is an AAN) in this case.

aus_economy %>% 
  model(
    mod = ETS(Pop)
  ) %>% 
  report()
Series: Pop 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.3266366 

  Initial states:
        l         b
 10.05414 0.2224818

  sigma^2:  0.0041

      AIC      AICc       BIC 
-76.98569 -75.83184 -66.68347 

A damped example: it’s leveling off.

An example with multiple series (and multiple models, best are selected automatically).

Lab session 14

Try foreasting Chinese GDP from global_economy using an ETS model.

Example: Australian holiday tourism. Observe all the unique combinations of models generated.

fit %>% pull(mod)
 [1] ETS(A,N,A) ETS(A,A,N) ETS(M,N,A) ETS(M,N,A) ETS(M,N,M) ETS(A,N,A) ETS(M,N,M)
 [8] ETS(M,N,A) ETS(A,N,A) ETS(A,N,N) ETS(M,N,N) ETS(M,N,M) ETS(A,A,N) ETS(A,N,A)
[15] ETS(M,N,A) ETS(A,N,N) ETS(M,N,M) ETS(M,N,A) ETS(M,N,M) ETS(M,N,M) ETS(M,N,A)
[22] ETS(M,N,N) ETS(M,N,A) ETS(M,N,M) ETS(A,N,A) ETS(M,N,A) ETS(M,N,A) ETS(M,N,A)
[29] ETS(M,N,A) ETS(A,N,A) ETS(M,N,A) ETS(M,N,A) ETS(A,N,A) ETS(A,N,N) ETS(A,N,A)
[36] ETS(M,N,M) ETS(M,N,A) ETS(M,N,M) ETS(M,A,A) ETS(M,A,A) ETS(M,N,A) ETS(M,N,A)
[43] ETS(M,N,A) ETS(A,N,A) ETS(M,N,A) ETS(A,N,A) ETS(A,N,N) ETS(M,N,A) ETS(A,N,A)
[50] ETS(M,A,A) ETS(M,N,M) ETS(A,N,N) ETS(M,N,M) ETS(M,N,M) ETS(A,N,A) ETS(M,N,A)
[57] ETS(A,N,A) ETS(A,N,A) ETS(M,N,A) ETS(M,N,M) ETS(M,N,A) ETS(M,N,M) ETS(A,N,A)
[64] ETS(M,N,A) ETS(M,N,M) ETS(M,N,N) ETS(M,N,A) ETS(M,N,A) ETS(M,N,A) ETS(A,A,A)
[71] ETS(M,N,A) ETS(M,A,A) ETS(M,N,M) ETS(M,N,M) ETS(A,A,N) ETS(M,N,A)

“Best” method chosen by AICc. See slide 39 of lecture 7. Paper for the method of “automatic forecasting” using AICc.

ETS won’t work well for data less than monthly level.

Lab session 15

Find an ETS model for the Gas data from aus_prouduction

fc %>% filter(.model == "log") %>% pull(.distribution)
 [1] t(N(5.3, 0.0033)) t(N(5.5, 0.0055)) t(N(5.5, 0.0081)) t(N(5.4, 0.0114)) t(N(5.3, 0.0158))
 [6] t(N(5.5, 0.0202)) t(N(5.6, 0.0254)) t(N(5.4, 0.0314)) t(N(5.3, 0.0389)) t(N(5.5, 0.0465))
[11] t(N(5.6, 0.0550)) t(N(5.4, 0.0645)) t(N(5.3, 0.0760)) t(N(5.6, 0.0876)) t(N(5.6, 0.1002))
[16] t(N(5.4, 0.1141))

Non-Gaussian forecast distributions

If you don’t want to believe normality. For instance, forecasts are too skewed (e.g. mean forecast not center of forecast distribution, or residuals not normally distributed).

For non-Gaussian resduals, you can bootstrap with generate(). Although this still assumed uncorrelated residuals. Can also do it within forecast() with bootstrap = TRUE. 500 samples is default number of samples, but you can changes this with times = N. See ?forecast.ETS().

Lab session

fit$ets
[1] ETS(A,A,A)

Bootstrapped forecast distributions

ETS handles trend, seasonality, but fails to handle cycles. For example:

# observe the model
pelt %>% 
  model(ets = ETS(Lynx)) %>% 
  pull(ets)
[1] ETS(A,N,N)
LS0tCnRpdGxlOiAiTGVjdHVyZSA3IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7cn0KbGlicmFyeShmcHAzKQpgYGAKCiMgRXhwb25lbnRpYWwgc21vb3RoaW5nCgpzZWUgcGRmIG5vdGVzIGZvciBkZXRhaWxzIG9uIGFkZGl0aXZlIGFuZCBtdWx0aXBsaWNhdGl2ZSBtb2RlbHMuCgpFVFMoZXJyb3IsIHRyZW5kLCBzZWFzb24pCgpFYWNoIGVycm9yLCB0cmVuZCBhbmQgc2Vhc29uIGNhbiBiZSBhZGRpdGl2ZS9tdWx0aXBsaWNhdGl2ZS9ub25lLCBvciBkYW1wZWQuCgpgYGB7cn0KYXVzX2Vjb25vbXkgPC0gZ2xvYmFsX2Vjb25vbXkgJT4lIAogIGZpbHRlcihDb2RlID09ICJBVVMiKSAlPiUgCiAgbXV0YXRlKFBvcCA9IFBvcHVsYXRpb24vMWU2KQoKZml0IDwtIGF1c19lY29ub215ICU+JSAKICBtb2RlbCgKICAgIEFBTiA9IEVUUyhQb3ApCiAgKQoKcmVwb3J0KGZpdCkKYGBgCgpFVFMoQSxBLE4pIDogQSBzdGFuZHMgZm9yICJhZGRpdGl2ZSIsIE0gZm9yICJtdWx0aXBsaWNhdGl2ZSIsIGFuZCBOIGZvciAibm9uZSIuIERhbXBlZCBpcyAiQWQiIG9yICJNZCIgZm9yIGFkZGl0aXZlIGFuZCBtdWx0aXBsaWNhdGl2ZSByZXNwZWN0aXZlbHkuIAoKJFxzaWdtYV4yJCBpcyB0aGUgdmFyaWFuY2Ugb2YgdGhlIHJlc2lkdWFscwoKYGBge3J9CmNvbXBvbmVudHMoZml0KQpjb21wb25lbnRzKGZpdCkgJT4lIGF1dG9wbG90KCkKYGBgCgpgYGB7cn0KZml0ICU+JSAKICBmb3JlY2FzdChoID0gMjApICU+JSAKICBhdXRvcGxvdChhdXNfZWNvbm9teSkKYGBgCgpEZWZhdWx0IGNob29zZXMgYmVzdCAod2hpY2ggaXMgYW4gQUFOKSBpbiB0aGlzIGNhc2UuCmBgYHtyfQphdXNfZWNvbm9teSAlPiUgCiAgbW9kZWwoCiAgICBtb2QgPSBFVFMoUG9wKQogICkgJT4lIAogIHJlcG9ydCgpCmBgYAoKQSBkYW1wZWQgZXhhbXBsZTogaXQncyBsZXZlbGluZyBvZmYuIApgYGB7cn0KZml0IDwtIGF1c19lY29ub215ICU+JSAKICBtb2RlbCgKICAgIG1vZCA9IEVUUyhQb3AgfiB0cmVuZCgiQWQiKSkKICApIApyZXBvcnQoZml0KQphdXRvcGxvdChmb3JlY2FzdChmaXQsIGggPSAyMCksIGF1c19lY29ub215KQpgYGAKCkFuIGV4YW1wbGUgd2l0aCBtdWx0aXBsZSBzZXJpZXMgKGFuZCBtdWx0aXBsZSBtb2RlbHMsIGJlc3QgYXJlIHNlbGVjdGVkIGF1dG9tYXRpY2FsbHkpLgpgYGB7cn0KZml0IDwtIGdsb2JhbF9lY29ub215ICU+JSAKICBtdXRhdGUocG9wID0gUG9wdWxhdGlvbi8xZTYpICU+JSAKICBtb2RlbChldHMgPSBFVFMocG9wKSkKCmZpdCAgCgpmb3JlY2FzdChmaXQsIGggPSA1KQpgYGAKCiMgTGFiIHNlc3Npb24gMTQKClRyeSBmb3JlYXN0aW5nIENoaW5lc2UgR0RQIGZyb20gZ2xvYmFsX2Vjb25vbXkgdXNpbmcgYW4gRVRTIG1vZGVsLgoKYGBge3J9CiMgY2FsY3VhbHRlIGdwIHBlciBjYXBpdGEgZm9yIGNoaW5hCmNoaW5hX2dkcCA8LSBnbG9iYWxfZWNvbm9teSAlPiUgCiAgZmlsdGVyKENvdW50cnkgPT0gIkNoaW5hIikgCgphdXRvcGxvdChjaGluYV9nZHAsIGdkcF9wYykKCiMgZml0IGV0cyBtb2RlbApjaGluYV9nZHAgJT4lIAogIG1vZGVsKAogICAgbW9kID0gRVRTKEdEUCkKICApICU+JSAKICBmb3JlY2FzdChoPTEwKSAlPiUgCiAgYXV0b3Bsb3QoY2hpbmFfZ2RwKQoKIyB0cnkgc29tZSB0cmFuc2Zvcm1hdGlvbnMKZml0IDwtIGNoaW5hX2dkcCAlPiUgCiAgbW9kZWwoCiAgICBldHMgPSBFVFMoR0RQKSwKICAgIGV0c19kYW1wZWQgPSBFVFMoR0RQIH50cmVuZCgiQWQiKSksCiAgICAjIGJveCBjb3ggdHJhbmZvcm1hdGlvbiB2aWEgR3VlcnJvIGlzIGFjdHVhbGx5IG5lZ2F0aXZlLCAKICAgICMgYnV0IG5ldmVyIHVzZSBzb2VtdGhpbmcgdGhhdCBsb3cgaW4gcHJhY3RpY2UKICAgIGV0c19iYyA9IEVUUyhib3hfY294KEdEUCwgMC4yKSksCiAgICBldHNfbG9nID0gRVRTKGxvZyhHRFApKQogICkgCgpmaXQgIAoKIyB3YXlzIHRvIGV4dHJhY3QgY29tcG9uZW50cwpyZXBvcnQoZml0KQpmaXQgJT4lIHNlbGVjdChldHMpICU+JSByZXBvcnQoKQpnbGFuY2UoZml0KQp0aWR5KGZpdCkKY29lZihmaXQpCgpmaXQgJT4lIAogIGZvcmVjYXN0KGg9IjIwIHllYXJzIikgJT4lIAogIGF1dG9wbG90KGNoaW5hX2dkcCwgbGV2ZWwgPSBOVUxMKQpgYGAKCkV4YW1wbGU6IEF1c3RyYWxpYW4gaG9saWRheSB0b3VyaXNtLiBPYnNlcnZlIGFsbCB0aGUgdW5pcXVlIGNvbWJpbmF0aW9ucyBvZiBtb2RlbHMgZ2VuZXJhdGVkLiAKYGBge3J9CmhvbGlkYXlzIDwtIHRvdXJpc20gJT4lIAogIGZpbHRlcihQdXJwb3NlID09ICJIb2xpZGF5IikKZml0IDwtIG1vZGVsKGhvbGlkYXlzLCBtb2QgPSBFVFMoVHJpcHMpKQpmaXQgJT4lIHB1bGwobW9kKQpgYGAKCiJCZXN0IiBtZXRob2QgY2hvc2VuIGJ5IEFJQ2MuIFNlZSBzbGlkZSAzOSBvZiBsZWN0dXJlIDcuIFBhcGVyIGZvciB0aGUgbWV0aG9kIG9mICJhdXRvbWF0aWMgZm9yZWNhc3RpbmciIHVzaW5nIEFJQ2MuCgoqKkVUUyB3b24ndCB3b3JrIHdlbGwgZm9yIGRhdGEgbGVzcyB0aGFuIG1vbnRobHkgbGV2ZWwqKi4gCgoKIyBMYWIgc2Vzc2lvbiAxNQoKRmluZCBhbiBFVFMgbW9kZWwgZm9yIHRoZSBHYXMgZGF0YSBmcm9tIGF1c19wcm91ZHVjdGlvbgoKYGBge3J9CmF1dG9wbG90KGF1c19wcm9kdWN0aW9uLCBHYXMpCgojIG5lZWRzIGEgdHJhbnNmb3JtYXRpb24KYXV0b3Bsb3QoYXVzX3Byb2R1Y3Rpb24sIGJveF9jb3goR2FzLCAwLjIpKQphdXRvcGxvdChhdXNfcHJvZHVjdGlvbiwgbG9nKEdhcykpCgpnYXMgPC0gYXVzX3Byb2R1Y3Rpb24gJT4lIAogIHNlbGVjdChHYXMpICU+JSAKICBmaWx0ZXIoUXVhcnRlciA8PSB5ZWFycXVhcnRlcigiMjAwNyBRNCIpKSAKYXV0b3Bsb3QoZ2FzKQoKIyBmaXQgbW9kZWxzCmZpdCA8LSBnYXMgJT4lIAogIG1vZGVsKAogICAgYXV0byA9IEVUUyhsb2coR2FzKSksCiAgICBkYW1wZWQgPSBFVFMoR2FzIH4gdHJlbmQoIkFkIikpLAogICAgbG9nID0gRVRTKGxvZyhHYXMpKQogICkKCmZpdCRhdXRvICAgICMgQUFBCmZpdCRkYW1wZWQgICMgTSBBZCBNCmZpdCRsb2cgICAgICMgQUFBCgpmYyA8LSBmaXQgJT4lIAogIGZvcmVjYXN0KGggPSAiNCB5ZWFycyIpCgpmYyAlPiUgCiAgYXV0b3Bsb3QoYXVzX3Byb2R1Y3Rpb24sIGxldmVsID0gTlVMTCkKCiMgem9vbSBpbgpmYyAlPiUgCiAgYXV0b3Bsb3QoCiAgICBmaWx0ZXIoCiAgICAgIGF1c19wcm9kdWN0aW9uLCAKICAgICAgUXVhcnRlciA+IHllYXJxdWFydGVyKCIyMDAwIFE0IikpLAogICAgbGV2ZWwgPSBOVUxMKQoKZmMgJT4lIAogIGFjY3VyYWN5KGF1c19wcm9kdWN0aW9uKQoKZ2xhbmNlKGZpdCkKYGBgCgpgYGB7cn0KIyBub3RlIHRoYXQgdGhlIGxvZyB0cmFuc2Zvcm1lZCBkaXN0cmlidXRpb25zIGFyZSB0cmFuc2Zvcm1lZCBub3JtYWwgZGlzdHJpYnV0aW9ucy4KIyBpdCBhcHBsaWVzIGEgYmlhcyBjb3JyZWN0aW9uIGZhY3RvciBpbiB0aGUgYmFja3RyYW5zZm9ybWF0aW9uIGF1dG9tYXRpY2FsbHkuID0pCmZjICU+JSBmaWx0ZXIoLm1vZGVsID09ICJsb2ciKSAlPiUgcHVsbCguZGlzdHJpYnV0aW9uKQpgYGAKCgojIE5vbi1HYXVzc2lhbiBmb3JlY2FzdCBkaXN0cmlidXRpb25zCgpJZiB5b3UgZG9uJ3Qgd2FudCB0byBiZWxpZXZlIG5vcm1hbGl0eS4gRm9yIGluc3RhbmNlLCBmb3JlY2FzdHMgYXJlIHRvbyBza2V3ZWQgKGUuZy4gbWVhbiBmb3JlY2FzdCBub3QgY2VudGVyIG9mIGZvcmVjYXN0IGRpc3RyaWJ1dGlvbiwgb3IgcmVzaWR1YWxzIG5vdCBub3JtYWxseSBkaXN0cmlidXRlZCkuCgpGb3Igbm9uLUdhdXNzaWFuIHJlc2R1YWxzLCB5b3UgY2FuIGJvb3RzdHJhcCB3aXRoIGBnZW5lcmF0ZSgpYC4gQWx0aG91Z2ggdGhpcyBzdGlsbCBhc3N1bWVkIHVuY29ycmVsYXRlZCByZXNpZHVhbHMuIENhbiBhbHNvIGRvIGl0IHdpdGhpbiBgZm9yZWNhc3QoKWAgd2l0aCBgYm9vdHN0cmFwID0gVFJVRWAuIDUwMCBzYW1wbGVzIGlzIGRlZmF1bHQgbnVtYmVyIG9mIHNhbXBsZXMsIGJ1dCB5b3UgY2FuIGNoYW5nZXMgdGhpcyB3aXRoIGB0aW1lcyA9IE5gLiBTZWUgYD9mb3JlY2FzdC5FVFMoKWAuCgoKIyBMYWIgc2Vzc2lvbiAKCmBgYHtyfQp2YyA8LSBhdXNfcmV0YWlsICU+JSAKICBmaWx0ZXIoU3RhdGUgPT0gIlZpY3RvcmlhIiwKICAgICAgICAgSW5kdXN0cnkgPT0gIkNhZmVzLCByZXN0YXVyYW50cyBhbmQgY2F0ZXJpbmcgc2VydmljZXMiKSAlPiUgCiAgc2VsZWN0KE1vbnRoLCBUdXJub3ZlcikKCnZjICU+JSBhdXRvcGxvdChUdXJub3ZlcikKCiMgYXBwbHkgdHJhbnNmb3JtYXRpb24gdG8gc3RhYmlsaXplIHZhcmlhbmNlIGFjcm9zcyB0aW1lIHNlcmllcwp2YyAlPiUgYXV0b3Bsb3QoYm94X2NveChUdXJub3ZlciwgMC4yKSkKCiMgZml0IG1vZGVsCmZpdCA8LSB2YyAlPiUgCiAgbW9kZWwoZXRzID0gRVRTKGJveF9jb3goVHVybm92ZXIsIDAuMikpKQoKZml0JGV0cyAKYGBgCgpCb290c3RyYXBwZWQgZm9yZWNhc3QgZGlzdHJpYnV0aW9ucwpgYGB7cn0Kc2ltIDwtIGZpdCAlPiUgZ2VuZXJhdGUoaCA9ICIzIHllYXJzIiwgdGltZXMgPSA1LCBib290c3RyYXAgPSBUUlVFKQpzaW0KCnZjICU+JSAKICBmaWx0ZXIoTW9udGggPj0gMjAwOCkgJT4lIAogIGdncGxvdChhZXMoeD1Nb250aCkpICsgCiAgZ2VvbV9saW5lKGFlcyh5ID0gVHVybm92ZXIpKSArCiAgZ2VvbV9saW5lKGRhdGEgPSBzaW0sIGFlcyh5ID0gLnNpbSwgY29sb3IgPSBhcy5mYWN0b3IoLnJlcCkpKSArCiAgZ3VpZGVzKGNvbCA9IEZBTFNFKSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSB5ZWFycXVhcnRlcihjKDIwMTAsIDIwMjMpKSwKICAgICAgICAgICAgICAgICAgeWxpbSA9IGMoMzAwLCA4MDApKQpgYGAKCgoKRVRTIGhhbmRsZXMgdHJlbmQsIHNlYXNvbmFsaXR5LCBidXQgZmFpbHMgdG8gaGFuZGxlIGN5Y2xlcy4gRm9yIGV4YW1wbGU6CmBgYHtyfQpwZWx0ICU+JSBhdXRvcGxvdChMeW54KQoKIyBoZXJlJ3MgYSBsb3VzeSBmb3JlY2FzdApwZWx0ICU+JSAKICBtb2RlbChldHMgPSBFVFMoTHlueCkpICU+JSAKICBmb3JlY2FzdChoPSIzMCB5ZWFycyIpICU+JSAKICBhdXRvcGxvdChwZWx0KQoKIyBvYnNlcnZlIHRoZSBtb2RlbDogZXJyb3IgaXMgQWRkaXRpdmUsIHRyZW5kIGlzIG5vbmUsIHNlYXNvbiBpcyBub25lLiAKIyBFVFMgbW9kZWwgZG9lc24ndCBrbm93IHdoYXQgdG8gZG8gd2l0aCBjeWNsaWMgZGF0YSB3aXRob3V0IHNlYXNvbmFsaXR5IQpwZWx0ICU+JSAKICBtb2RlbChldHMgPSBFVFMoTHlueCkpICU+JSAKICBwdWxsKGV0cykKYGBgCgo=